% This file take the steady-state the file "todynaremone.mat"
% It writes the dynare file "code_dynare_rules.mod" which is the dynamic model
% to use perturbation techniques, using fiscal rules
% it doubles check that the steady state is OK, thanks to the "resid" function

% 


clear
warning('off','all');
isOctave = exist('OCTAVE_VERSION', 'builtin') ~= 0;

isOctave && warning('off', 'Octave:array-as-logical');

load ../toDynare

tauk = eco.tauk;


alpha = eco.alpha(1);
beta = eco.beta;
delta = eco.delta(1);
kappa = eco.kappa;
chi = eco.chi;
phi = eco.phi;
G = solution.G;
coefB = 0.1;


rho_u = 0.95;
G0 = .01*(1-rho_u); 


taut =  (( 1.0/eco.phi + 1.0)*eco.tau)/( 1.0/eco.phi +1.0 - eco.tau);

file = 'code_dynare_RA_rules.mod';
fid  = fopen(file, 'w');
formatSpec = '%25.20f';

GsY = solution.G/solution.Y;
chi = 1;

BsY = solution.B/solution.Y;

dif = 1;


r = 1/beta - 1;
FK = r/(1-tauk);
KsL = ((FK+delta)/alpha)^(1/(alpha - 1));
FL = (1-alpha)*KsL^alpha;
w = kappa*FL;
l =  (chi*taut*w)^( (1/phi+1 +taut)/((1+1/phi)*(1+1/phi)) );
K = KsL*l;
Y = K^alpha*l^(1-alpha);
I = delta*K;
c = Y*(1-GsY) -I; 
G = GsY*Y;
x = c - (1/(chi*(1+1/phi)))*l^(1+1/phi);
a = (x - (1/(chi*taut))*l^(1+1/phi))/r;
B= a-K;


% using chi to match B/Y
while abs(dif)>10^-8

r = 1/beta - 1;
FK = r/(1-tauk);
KsL = ((FK+delta)/alpha)^(1/(alpha - 1));
FL = (1-alpha)*KsL^alpha;
w = kappa*FL;
l =  (chi*taut*w)^( (1/phi+1 +taut)/((1+1/phi)*(1+1/phi)) );
K = KsL*l;
Y = K^alpha*l^(1-alpha);
I = delta*K;
c = Y*(1-GsY) -I; 
G = GsY*Y;
x = c - (1/(chi*(1+1/phi)))*l^(1+1/phi);
a = (x - (1/(chi*taut))*(l^(1+1/phi)))/r;
B= a-K;

dif = BsY - B/Y;
chi = chi + 0.001*dif;
end


cc = x - ( (1/(chi*taut))*l^(1+1/phi)  + r*a ) 

%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%     variables and params     %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

str = ['var\n'] ; fprintf(fid, str);

str = ['u r w FK FL B K TFP Ltot Ctot Y I It\n'] ; fprintf(fid, str);
str = ['ut Kt Ct Bt wt rt taut tauk Yt Lt kappa G l T a x ;\n'] ; fprintf(fid, str);


str = ['varexo eps; \n\n'] ; fprintf(fid, str);
str = ['parameters\n'] ; fprintf(fid, str);
str = ['beta alpha phi abar delta gamma chi rho_u Gss;\n\n'] ; fprintf(fid, str);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%        initialisation        %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


    
    
str = ['alpha','   = ',num2str(alpha,formatSpec),';\n']; fprintf(fid, str);
%str = ['beta','    = ',num2str(eco.beta,formatSpec),';\n']; fprintf(fid, str);
str = ['beta','    = ',num2str(beta,formatSpec),';\n']; fprintf(fid, str);
str = ['phi','     = ',num2str(eco.phi,formatSpec),';\n']; fprintf(fid, str);
str = ['abar','    = ',num2str(eco.a_min,formatSpec),';\n']; fprintf(fid, str);
str = ['delta','   = ',num2str(eco.delta,formatSpec),';\n']; fprintf(fid, str);
str = ['gamma','   = ',num2str(eco.sigma,formatSpec),';\n']; fprintf(fid, str);
str = ['rho_u','   = ',num2str(rho_u,formatSpec),';\n']; fprintf(fid, str); %HERE
str = ['Gss','     = ',num2str(G,formatSpec),';\n']; fprintf(fid, str);
str = ['chi','     = ',num2str(chi,formatSpec),';\n']; fprintf(fid, str);

%%
equa = 1; % Numbering equations
tol = 0*10^-8;  % threshold for transition, to avoid to consider very small transitions

str = ['\n\n'] ; fprintf(fid, str);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%       model definition       %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


str = ['model;\n\n']; fprintf(fid, str);

    
str = ['x = (1/(chi*taut))*l^(1+1/phi)  + T + (1+r)*a(-1)-a; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['(x^-gamma) =  beta*(1+r(+1))*(x(+1)^-gamma); %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
%%%%%%%% Resource contraint
str = ['G = Gss*(1+u); %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['I = K - (1-delta)*K(-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['G  + Ctot + I = Y;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

str = ['a = B + K; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;


%%%%%%%% Production side
str = ['FL = (1 - alpha)*TFP*(K(-1)/Ltot)^alpha; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['FK = alpha*TFP*(K(-1)/Ltot)^(alpha-1)-delta; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['u = rho_u*u(-1) + eps;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['TFP = 1; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Y = TFP*K(-1)^alpha*Ltot^(1-alpha);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

%%%%%%%% Wage
str = ['l =  (chi*taut*w)^( (1/phi+1 +taut)/((1+1/phi)*(1+1/phi)) ); %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;



%%%%%%%% Aggregate labor supply
str = ['Ltot = l; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

%%%%%%%% Aggregate consumption
str = ['Ctot = (x +(( l^(1 + 1/phi)) /(chi*(1/phi+1))) );  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;


%%%%%%%% Taxes
str = ['tauk = 1-r/FK;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['kappa = w/FL;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;


%%%%%%%%% Fiscal system

str = ['tauk  = steady_state(tauk);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['kappa = steady_state(kappa);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['taut = ',num2str(taut,formatSpec),';  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;


str = ['T = -',num2str(coefB,formatSpec),'*(B - steady_state(B)) ;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;



%%%%%%%%%%%%%%%% log-deviation variables %%%%%%%%%%%%%%%%
% There are denoted with 't': for variable x, xt = (x -x_ss)/x_ss, where x_ss is the steady state value of x.

str = ['ut = 100*u;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Kt = 100*(K/',num2str(K,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Ct = 100*(Ctot/',num2str(c,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Bt = 100*(B/',num2str(B,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['wt = 100*(w/',num2str(w,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['rt = 100*(r - (1/beta - 1));  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Lt = 100*(Ltot/',num2str(l,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Yt = 100*(Y/',num2str(Y,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['It = 100*(I/',num2str(delta*K,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['end;\n\n'] ; fprintf(fid, str);

str = ['%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n\n'] ; fprintf(fid, str);
%%%%%%%%%%%%%%%% end of the model part


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%     Steady-state model       %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



str = ['initval;\n\n'] ; fprintf(fid, str);
    
str = ['l     = ', num2str(l,formatSpec),';\n'] ; fprintf(fid, str);
str = ['r     = 1/beta - 1;\n'] ; fprintf(fid, str);
str = ['w     = ', num2str(w,formatSpec),';\n'] ; fprintf(fid, str);
str = ['FK    = ', num2str(FK,formatSpec),';\n'] ; fprintf(fid, str);
str = ['FL    = ', num2str(FL,formatSpec),';\n'] ; fprintf(fid, str);
str = ['K     = ', num2str(K,formatSpec),';\n'] ; fprintf(fid, str);
str = ['x     = ', num2str(x,formatSpec),';\n'] ; fprintf(fid, str);
str = ['Ltot  = ', num2str(l,formatSpec),';\n'] ; fprintf(fid, str);
str = ['TFP   = ', num2str(1),';\n'] ; fprintf(fid, str);
str = ['B     = ', num2str(B,formatSpec),';\n'] ; fprintf(fid, str);
str = ['a     = ', num2str(a,formatSpec),';\n'] ; fprintf(fid, str);
str = ['u     = ', num2str(0),';\n'] ; fprintf(fid, str);
str = ['tauk  = ', num2str(tauk,formatSpec),';\n'] ; fprintf(fid, str);
str = ['taut   = ', num2str(taut,formatSpec),';\n'] ; fprintf(fid, str);
str = ['Ctot  = ', num2str(c,formatSpec),';\n'] ; fprintf(fid, str);
str = ['Y     = ', num2str(Y,formatSpec),';\n'] ; fprintf(fid, str);
str = ['kappa = ', num2str(kappa,formatSpec),';\n'] ; fprintf(fid, str);
str = ['I  = delta*K;\n'] ; fprintf(fid, str);

str = ['ut    = 0;\n'] ; fprintf(fid, str);
str = ['Kt    = 0;\n'] ; fprintf(fid, str);
str = ['Ct    = 0;\n'] ; fprintf(fid, str);
str = ['Bt    = 0;\n'] ; fprintf(fid, str);
str = ['wt    = 0;\n'] ; fprintf(fid, str);
str = ['rt    = 0;\n'] ; fprintf(fid, str);
str = ['Lt    = 0;\n'] ; fprintf(fid, str);
str = ['Yt    = 0;\n'] ; fprintf(fid, str);
str = ['It    = 0;\n'] ; fprintf(fid, str);
str = ['T  = 0;\n'] ; fprintf(fid, str);


str = ['G     = ', num2str(G,formatSpec),';\n'] ; fprintf(fid, str);

str = ['end;\n\n'] ; fprintf(fid, str);
%%%%%%%%%%% end of the steady state part

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%   Steady-state computation    %%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

str = ['resid;\n\n'] ; fprintf(fid, str);
str = ['steady(nocheck);\n\n'] ; fprintf(fid, str);
 
%str = ['check;\n\n'] ; fprintf(fid, str);



    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%       Aggregate shock        %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 str = ['shocks;\n var eps =',num2str(G0,formatSpec),';\nend;\n'] ; fprintf(fid, str);



    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%         Stoch simul          %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
%~ str = ['stoch_simul (order=1,irf=200) ut Bt tau kappa tauk Yt Ct Kt Lt TLinc TKinc G Y BoYt mut Epsi corre;'] ; fprintf(fid, str);
%~ str = ['stoch_simul (order=1,irf=200) ut Bt tau kappa tauk Yt Ct Kt Lt G Y BoYt mut;'] ; fprintf(fid, str);


str = ['stoch_simul (order=1,irf=500) ut Bt taut kappa tauk Yt Ct It Lt;'] ; fprintf(fid, str);

isOctave && fflush(fid);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%        Launch Dynare         %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    command = ['dynare ',file, ' noclearall'];
    eval(command);

%name = ['fig',num2str(rho_u),'.mat'];
%save (name,'ut_eps','Bt_eps','tau_eps','kappa_eps','tk_eps','Yt_eps','Ct_eps','Kt_eps','Lt_eps','TLinc_eps','TKinc_eps','G_eps','Y_eps','BoYt_eps','mut_eps','oo_')


save ('fig_RA95_C4.mat','ut_eps','Bt_eps','Yt_eps','It_eps','Lt_eps','Ct_eps')




